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QUEST: A Quasi-Equilibrium Transport 


Code for High-Beta Systems 


I. INTRODUCTION 


There are an increasing number of high-beta systems where the plasma dynamical time, 
defined as the time it takes the plasma to adjust to a force imbalance, is much shorter than the 
plasma lifetime. In reversed-field configurations, pinches and, perhaps, mirrors, the 
confinement time will be determined by diffusion processes, either classical or anomalous. Also, 
in liner experiments there is compressional heating on the slow (liner-speed) time scale. 
Diffusion and liner times are much longer than plasma dynamical times, which are typically the 
sonic and Alfven transit times. Thus to avoid prohibitively stringent time-step limitations it is 


necessary to develop a code which operates on the slower time scales. 


It is well known that in low-8 systems such as tokamaks, one can eliminate the dynamical 
time scale by neglecting the diffusion of the toroidal magnetic flux to lowest order and thus 
obtain an algebraic expression for the diffusion velocity and coefficient. However, in high-B 
(B = 1) systems, the particle and magnetic flux diffusion are comparable and it is not possible 
to write an expresssion for the diffusion velocity and coefficients in closed from (see, e.g. 


Liewer and Davidson!!)), 


But the diffusion processes are slow compared with the time for the plasma to reach force 
balance, so that the plasma essentially remains in equilibrium throughout the process. In a 
manner similar to that employed earlier to study impurity transport,!2) the QUEST code makes 
use of this to eliminate the dynamical time scale. Once diffusion has taken place so that the in- 


itial equilibrium plasma configuration moves slightly out of force balance, the code returns the 
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configuration to equilibrium by assuming the plasma relaxes adiabatically and conservatively, 
i.e., with the particle number, magnetic flux and entropy held constant. The plasma is then al- 
lowed to diffuse again, and the relaxation repeated. This eliminates the need to solve the 
momentum equation, and thus eliminates the dynamical time scale. Thus it becomes feasible 
to model the evolution of the system over times as long as 10°—104 sound or Alfven transit 


times, which would require as many as 10° time steps with conventional algorithm. 


Compressional heating is included in the same way. After the plasma is compressed 
slightly by, say, inward motion of the outer wall, the plasma is again slightly out of force bal- 


ance. Again, the configuration is returned to equilibrium assuming an adiabatic relaxation. 


Thus far the code has been used to study the classical diffusion in the SEEBIE reversed 
magnetic field experiment at NRL." This experiment has comparable axial and azimuthal mag- 
netic fields. The code has also calculated successive states of an ideal lossless (i.e., with all tran- 


sport processes turned off) reversed field configuration. 


In the future, the code will be used to study strong compressional heating of such 
configurations both for the present parameters and for reactor parameters. Anomalous tran- 
sport coefficients will be added. It is also planned to study the stability of such configurations 
by using the equilibrium configurations generated by QUEST as input to MHD and tearing 


mode stability codes.'*! (Thus QUEST can stand for quasi-equilibrium stability and transport.) 


The organization of this report is as follows. In Section II, the equations solved by the 
code are presented, and the method of solution summarized. In Section IJ], we describe how to 
run the code. The input parameters and interpretation of the output are discussed. In Section 
IV, a description of each subroutine, its purpose and the solution methods used are described. 
Section V briefly describes the results obtained when measured SEEBIE profiles are used to ini- 


tialize studies of adiabatic compression and classical transport separately. Appendix A contains 


—— 
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a description of the iteration method used to relax the grid to equilibrium. The code is listed in 


Appendix B. Sample input and output appears in Appendix C. 


II. SUMMARY OF CODE 


In this section, the equations solved by QUEST are presented and the method of solution 


described briefly. 


A. Basic Equations 


The evolution of a quasi-neutral plasma of two species [electrons (e) and ions (/)] is fully 


described by the equations 


n=n,=nN,, (1) 
nm Sy: Gee Le KB~ Up. (2) 
i or c ? 

V x B=4nJ/c; (3) 
oe cV x E; (4) 

or 
P= nk p( T, + T); (5) 

3 98 3 

a (nT) + V -(nT,v) + 17,97 -v=Q, + (6) 


V -(KEViT. + KE VT), a =i, e; 

E+ Lyx B=nJ,+ ni, (7) 
where we have assumed a scalar pressure p, and Q, is a heating term which includes resistive 
heating. The units are cgs through out this report and in the code except for the temperatures, 
which are in electron volts. Boltzmann’s constant is kg = 1.6 x 10~!? ergs/eV. The subscripts | 
and || refer to directions perpendicular and parallel to the magnetic field B. The atomic number 


Z of the ions is taken to be unity. 


We assume a cylindrically symmetric system which is infinite in the axial (z) direction. 


The equations then become, using Braginskii's'*! transport coefficients 
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1 
= ke (nid? + m2) — OQ, + Qen sas 
= nae [D2J? + D3J§ + 2DJoJ,] — Q, + Qea- 
B 


In obtaining the expression for the ion heat conductivity in Eq. (22), we have assumed 
the ions are magnetized, w,7; >> 1. The Q,, are any additional heating or cooling terms 


which may be included. 


Equations (7-25) are the complete set of equations solved by the code. They are not, 
however, solved by a time integration, but rather in the quasiequilibrium approximation as 


described below. 


In the absence of any diffusive processes, the right hand sides of Eqs. (6-7) vanish and 


Eq. (8) and Eqs. (10-12) can be written in the form 


Sry tk & 
Vlas a ar 


91 (rv,F) =0 (26) 


where 


F, =n [Eq. (8)]; 


Fy=n8T —— {Bq.(10)); 
F,= Br = Wg ADD): 


The four quantities F, — F4 times the volume element 27rdr are the usual conserved quanti- 
ties: particle number, entropy, and 6- and z-projections of magnetic flux, respectively. Conser- 
vation of energy can be derived analytically from Eqs. (8-12) if the transport terms are 
dropped, but, in general, need not hold for numerical approximations that conserve the four 


quantities associated with F,; — Fy. 


B. Method of Solution 


In the code Eqs. (7-25) are not integrated on the dynamical time scale. Instead, as dis- 


cussed in the introduction, it is assumed that the processes which destroy the equilibrium [the 
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right hand sides of Eqs. (8-12)] operate on a time scale long compared to the time to bring the 
plasma to force balance. This allows Eqs. (7-25) to be solved in two steps. The plasma starts 


in equilibrium 


4 Be - LB) = Se (27) 
Thus, from Eq. (9) v, = 0. Equations (7-25), with v, = 0 (i.e., only the right hand sides of 
Eqs. (8-12) contribute to the time rates of change of the dependent variables) are then 
advanced using a fully time- and space-centered implicit scheme. The diffusion results in a 
force imbalance across the grid (J x B = Vp). The configuration is brought back to equili- 
brium by an adiabatic relaxation of the grid: each cell has associated with it a value of the den- 
sity n,, an entropy function 7,,/n//°, and values of B,,/r, and B.,. Each of these four quanti- 
ties multiplied by the volume, dV, = 2mr,Ar, where Ar, is the width of j cell, is held constant 


while an iteration on the cell widths Ar, is performed until pressure balance is again restored. 


(The numerical procedure for doing this is described in Appendix A.) This is equivalent to 


finding the steady state solution 


Ov, 
ae Vv, = | to Eq. (27) plus the momentum equation (9). 


In addition, the code allows a third step, a compressional heating. The outer radius R of 


the cylinder moves according to 


aR 


Mi (28) 
Thus the outer gridpoint is moved in a time step, Av, by an amount 
AR = Vivant (29) 


This also results in a force imbalance across the grid. The grid is again brought back to equili- 


brium by mean of the adiabatic relaxation described above. 


The boundary conditions currently in the code specify 7, = 7, = 9 at the axis and 


at the wall, so that the material pressure vanishes at both ends of the coordinate system. The 


oF 
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wall is taken to be a perfect thermal and electrical conductor and impervious to material flow, 
so that v, = V,,,, there. The boundary conditions are used in initializing the code, solving the 
transport and pressure equilibration equations, and outputting the data. If different conditions 


are implemented in their place, changes should be made consistently throughout the code. 


III. RUNNING THE CODE 


In this Section, we discuss how the code is run. We describe first the input parameters, 
then the output. The code itself is treated here as a black box. Units are cgs throughout except 


for temperatures, which are in eV. The ion species is assumed to be deuterium. 


A. Input data 


Input is entered through two namelists: INPUT, which contains the parameters for the 
equilibrium conditions, and KNOB, which contains logical switches for various options in the 


code. 


1. Namelist INPUT 


GAM — the adiabatic constant y. Currently this must be set to 5/3. 
DRG — original grid spacing 
ETA — parameter specifying initial ion and electron temperature 


ratio, n = 7,/T, 


BZO — initial axial field at the wall 
RRIGHT — _ initial position of the wall 
RWALL — _ initial wall radius (should equal RRIGHT) 
RMIN — radius of wall at peak compression 
TMID — time of peak compression 
TO — initial electron temperature 
7 
= Ee ee ace ee Tene % = - pert WS Srp aeees a ne 
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DT — time step (Since code is implicit, there is no stability 
requirement on DT. The user should select it for a good 
representation of the diffusional and compressional effects.) 
BTO — parameter determining initial A, field at the wall (see Section IV B). 


NSTOP — maximum number of time steps 


NPRINT — print output and page plots every NPRINT steps f 
NPLOT -—_ produce calcomp plots every NPLOT steps if NPLOT is a 
multiple of NPRiNT 


2. Namelist KNOB 


LDSP — if true, dissipative processes are allowed (i.e., right hand sides 
of Eqs. (7-25) are firite); otherwise, only adiabatic compression 


LWALL — _ if true, the outer wall is moved in according to 


d RRIGHT — RMIN .. 
at R yait Vwalt TMID ; if false 


wall remains stationary. 
LPLOT -— _ if true, generates page plots and calcomp plots 
LBOOK -—_ currently must be true; eventually will allow a reverse « 
bias theta-pinch option 
LTURB -—_ currently must be false; eventually will allow for 
turbulent transport 
The code is now set to simulate the experiment SEEBIE and the program set for 28 grid points. 


(This can be changed by a macro instruction.) 
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B. Output 

Various grid quantities are printed out: R, the location of the grid points; B, the axial 
magnetic field (B,); N, the particle density; TE, the electron temperature; TI, the ion tempera- 
ture; PTOTAL = nk,(T. + 7,) ; BP, the magnetic field pressure pg, defined as 


B? o _, Bet) 
fais + far ————* 


a 8 A 4rr 


K, the total plasma plus magnetic field pressure in each cell (this must be the same in each cell 
or the grid relaxation has failed); BT, the azimuthal field B,; JZ, the axial current density; JT, 
the azimuthal current density, DELVZ, defined as —J./ne, DELV, defined as —J/,/ne, RESIST, 


the perpendicular resistivity 7). 


Also printed out are various extensive parameters characterizing the entire configuration: 


NTOTAL —_ total number of particles 
TOTE — total (plasma plus magnetic) energy 
BENER -— _ total magnetic energy 
PENER- — _ total plasma energy 
ITOTZ — total axial current 
ITOTT — total azimuthal current 
© BZTOTT — _ total axial flux 
BITOT -—_ total azimuthal flux per cm 
DIAM — diamagnetic signal ® measured at r = R(23) (set for SEEBIE), 
R (23) 
®=27 f r dr BAr) 
0 

BTRAT —_ B,(23)/B2(23) where B? is the value at ¢ = 0. 
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IV. DESCRIPTION OF CODE AND NUMERICAL METHODS USED 


This section contains a description of each subroutine in the code. The purpose of each 
one is presented and the numerical methods used are described. The code itself is listed in 


Appendix B. 


A. LINER — the main program ¢ 


E: Initializes various constants such as species charge and mass. (Throughout code, 
species I is electrons and 2 is deuterium ions.) 

z Reads the input parameters through the namelists described in Section III A. 

3; Calls the subroutine INCON which initializes the grid quantities and sets up an 
equilibrium configuration. 

4. Calls the subroutines FLUX and LAGRAG to make sure the initia! equilibrium 
configuration is a "finite-differenced” equilibrium. 


i Calls the subroutine STEP which now takes over control of the time evolution of the 


initial configuration. 


B. Subroutine INCON . 
Function: initializing the grid quantities. 
There are NRP grid points and NR = NRP-1 cells. Quantities can be defined on the grid 
point fe.g., RG(W) = 1/2 IRV) + RG + 1). Quantities ending in G (RG, DENG, etc.) are 
defined on the grid. Otherwise they are cell quantities (R, DEN, etc.). Most variable names are 
self explanatory. DR(J) is the array of cell sizes, DRU) = RGU + 1) - RG(). The initial cell 
quantities are stored in arrays with the same names but with the suffix OLD, e.g., BLOLDU) 
= BZ(J) at 1 = 0. 
10 
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In this version, first the B, and B. fields are calculated: 


B,(r) = BTO dr reo’ ~ RIP)2/DRIZ 
r 


BLO + BOL). + BZOL) y~(r ~ RIP)2/DRIT 
-DRIT Sy ar « 


where BTO, BZO, DRJZ, DRJT and RJP are input parameters. BZOL is calculated in INCON 


B.(r) = — BZOL + ——_——— 


: so that the plasma pressure is zero on axis. BZO is approximately the BZ field at the wall, RJP 
is the radius at which the current peaks, DRJZ and DRJT are related to the widths of the 


current profiles and BTO is related to the B, field at the wall. 
Once the magnetic field profiles are calculated, the plasma pressure profile is calculated 
from force balance: 


B} + B2 t de Be(r) 
oO im oe : +f dr 4 


Or if 8a Pa 4n P 


The temperature profiles are specified and the density profiles then calculated from the pressure 


and temperature profiles. 
C. Subroutine STEP 


Function: controls the time loop. Successive actions are as follows: 


1. Update NSTEP, the number of time steps taken and CLOCK, the time, and calculate 


the adiabatic invariants by calling the subroutine FLUX. 


2. Move the outer wall in (if LWALL is true; see Section IV A) by calling the subrou- 
tine WALL. 

3. Readjust the grid adiabatically to achieve force balance by calling the subroutine 
LAGRAG. 

4. Calculate the new grid quantities from the adiabatic invariants and the new cell sizes. 


This is done in by calling the subroutine OUTPUT. 


tant ry 
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5. Allow diffusion and dissipation if LDSP is true and then returns the grid to force bal- 


ance by calling the subroutine FLOW. 


6. Produce output if time step condition is met by calling OUTPUT. 


These six steps are then repeated NSTEP times. 


¢ 
D. Subroutine FLUX 
Function: calculate the four adiabatic invariants [see Eq. (26)] which are held constant dur- 
ing the grid relaxation to achieve pressure balance. The four quantities actually calculated 
are 
MASS (5) — MASS in the Jrhcell =m.n (J) x m(RGU + 1)? — (RGGQ))4 
§,() — entropy function, $,J) = P,(J)/[m,n())*? 
BZFLX (J) — axial flux = BZ(J) x m [((RG(W + 1 ))? — (RG(J))7) 
BTFLX (J) -- azimuthal flux BT() x 7[RGQU + 1) — RG(J)] 
E. Subroutine WALL 
° 
Function: computes the wall radius RWALL as a function of time (denoted CLOCK) 
according to RWALL = RMIN + (RRIGHT — RMIN) x ABS (CLOCK — TMID) / 
TMID. The parameters are described in Section 'V A. 
F. Subroutine LAGRAG 
Function: relax the grid to force balance by using an iteration process developed by 
Boris,'2! described in Appendix A. It calls the subroutine DDELR. 
12 
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G. Subroutine FLOW 
Function: Sets up the tridiagonal matrices to solve the equations for advancing 

T., T,, Bo, B, (eqs. (10-12) with v, = 0)]. These equations are solved implicity by inverting 
the tridiagonal riiatrices. Furthermore, since the transport coefficients (K,, 7, etc.) depend on 
T,, and B, an iteration is performed on B, 7, and the transport coefficients until the desired 
accuracy is achieved. This is done with the following procedure at each time step: 
1. Original values of B, 7, and the transport coefficients are stored (these matrices have 

names prefixed with an O, e.g., OBZ (J) and OKAPE(J). 
2: An iteration loop is preformed to achieve the desired accuracy. 

a. The subroutine COLL is called to calculate the transport coefficient using the current 


(iterating) values of B and 7, 
b. The tridiagonal matrtices are evaluated using transport coefficients of the form 
K.= ; (K" + K2) 
where K{ is the current transport coefficient [(calculated in step (2a)] and K2 is the-origi- 


nal, stored value [(calculated in step (1)]. 


c. The new values of 7, and B are calculated by inverting the tridiagonal matrices (by cal- 


ling the subroutine TRID). 


d. The new temperature is compared to the previous value to see if the desired accuracy 


has been achieved. 


H. Subroutine Coll 


Function: Calculates the classical transport terms appearing in eqs. (10-12): 


electron-ion collision frequency 
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TAU (1, 1) = 1,, defined in Eq. (18) 
ion-ion collision frequency 
TAU (1,2) = 1,, defined in Eq. (19) 
Perpendicular and parallel resistivities 


RESPRP (1) = »,, defined in Eq. (16) 


RESPAR(I) = 7», defined in Eqn. (17) . 
Magnetic field diffusion coefficient 

BDONE(I) = Dj, defined in Eq. (13) 

BD2(1) = D,, defined in Eq. (14) 

BD3(1) = D3, defined in Eq. (15) 
ion heating rate 

HT (1, 2) = Q, (with Q,, = 0), defined in Eq. (24) 
electron heating rate 

HT(1, 1) = Q, defined in Eq. (25) [In the Seebie run, there was a Q,, which modeled 

thermal conduction losses out the end] 
cross field thermal conductivity 

’ 


KAPE(I) = K,, defined in Eq. (21) 


KAPI(I) = K,, defined in Eq. (22) 


Vv. EXAMPLES 


A, Adiabatic Compression 


We can turn off the transport effects and heat sources and sinks by setting 


Q, = K,=7 = 0. The total energy is given by 


2 
Woy = 2 fo" rd nkglT, + 7) + P|. (30) 
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The rate of change of W,,, is now just equal to the work done on the system by the imploding 
liner, 


Woy = AR yu MV wai? (31) 


where the total pressure Pis given by 

P = nk,(T, + T,) + B?/8n, (32) 
“where the first term should be small or zero. The total pressure clearly increases with decreas- 
ing R,. There is an approximate power-law dependence which becomes exact in the trivial 


cases of pure magnetic flux or pure plasma'®!, 


For initial conditions, we chose the equilibrium profiles shown in Fig. 1. They were 
found by numerically fitting the profiles in the model used by Sethian, et. al.,'7! which were 
obtained in turn by an analytic fit to observations. There is, of course, no compression in this 
experiment, and the parameters of the observed configuration are unsuitable for a slow liner 
design (temperatures are of order 5-10 eV). For a test of the code and a preliminary exercise 
in preparation for realistic calculations using electron-beam-induced reversed-field 
configurations with scaled-up parameters, however, these initial conditions are a reasonable 


choice. 


Using V4, = 10* cm/sec in Eq. (28), we find the profiles shown in Fig. 2 after 0.3 msec 
and in Fig. 3 after 0.6 msec (i.e., after radial compression ratios of 1.75 and 7.0, respectively). 
The latter state has a maximum total pressure 915 times greater than in the initial state, and 
maximum B, and T of 30.4 and 16.6 times those at ¢ = 0. Note that the profiles have 
steepened considerably. This trend is expected to reduce the margin of MHD stability origi- 
nally present, and (for sufficiently high compression) to destabilize the equilibrium. Even 
before this effect is noticeable, the steeper gradients are likely to drive microinstabilities which 
enhance particle and heat diffusion, thus tending to reduce the effectiveness of the compression 


in heating the plasma to thermonuclear conditions. 


15 
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B. The SEEBIE Experiment 

Returning to the basic system of equations (1-7), we adopt the Braginskii (classical) tran- 
sport coefficients, Eqs. (13-25), and attempt to study the evolution of the initial profiles shown 
in Fig. 1, in the absence of compression. Two important competing effects are thereby intro- 
duced: generation of heat from decay of the magnetic fields (ohmic heating), and cross-field 
thermal transport, primarily carried by the ions. 

In order to obtain agreement with the observed"! decay of this configuration, it is neces- 
sary to include a model of axial heat losses by electron conduction. Though, strictiy speaking, 
this requires a two-dimensional model, a crude representation is obtained by adding to Eq. (25) 


a term of the form 


Qe ~ —2nT,¥,/I, (33) 


where v, = (kgT,/m,)'/? is the electron thermal speed, and / = 40 cm. The resulting values of 
the diamagnetic signal ® and of By are compared with those measured by probes at 
approximately 6 cm radius. Although the oscillations observed experimentally in AB. are not 
reproduced, the overall agreement is quite good, as is to be expected from a fluid description of 
this stable, cold plasma. | 

Appendix C contains the initial code parameters used to perform this calculations and that 
described in Section V A above, as well as the printed output (described in Section III B) ini- 


tially and at the end of each run. 
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Appendix A 


DETAILS OF THE PRESSURE BALANCE ALGORITHM 


Here we show how the pressure balance Eq. (27) can be used to determine the motion of 


the numerical Lagrangian cells. We first take the reference position ro equal to zero and write 


Eq. (27) as 
Br)  ¢, Be _ 
p(r) + = + Sa a = K, (A-1) 
where 
2 
K = p(0) + BO) 
87 


We write (A-1) in finite difference form 


Br, 41/2) 1 J BE ry 1) a Fis 
+ ee a + — —. 
PArisii2) 80 an S (rg a) + r4)/2 
] BA (ri 41) (Cr 41 <= r,)/2 
tn ees. = K. (A-2) 


Pressure balance adjustment occurs under the condition that the magnetic flux, number of par- 
ticles, and entropy, within each cell, are conserved. We now estimate p and Bin (A-2) by refer- 


ence to the conserved quantities. The magnetic fluxes are related to the fields by 


(5@),, / 
Bylr,44/) = i (A-3) 
(410°C 
(5®.) 41/2 
BAr,4\;.) = —— .. (A-4) 
ee OR = 
To eliminate p from Eq. (A-1) we note that 
p= TS." (A-5) 


SS 
Heenan 


—mesmnst 
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where a refers io the species for which the entropy function S has been defined. When only 
One ion species is present (as in the versions of the code described here), we write p, = nm, 
and p, = nm,.. For multispecies problems we introduce the number density fraction of species 


a, f= n,/n, whence 


Pa = PLM, 2 fama). (A-6) 
Now the total density in cell / + 1/2 is 
(8M) .4/2 
(aaa) = a, (A-7 
PX 41/2 eh, = 5 ) 


where (8M) ,,/) is the total mass in the cell. Substitution of (A-7) into (A-6) and the result 
into (A-5) allows us to write (A-2) in terms of the cell boundary positions and quantities held 


constant in the pressure balance adjustment: 


(8M) is2 fate |" (8) 2) (8@.) 3.4/2 
D Sa(ri41/2) So eee jT<-., , + aE A + Pine — 
a eras — Fi) (2 fom) Sarria, — rp)? Barr Ay — 7/7) 
; 
‘=! (8@,) ? (5@,) 2.1, 
ia eee L C a8 ak. (A-8) 
2m gel lest — Fé ae th) <7; 


Consider now a system for which it is known that a pressure balance existed a short time previ- 
ously. If, in the subsequent time interval the "constants" of the motion changed because of dis- 
sipative effects, then a new pressure balance with the "new" constants of the motion is required. 
Even neglecting dissipative effects a new pressure balance is required when a boundary moves 


or flux is added to the system. 


The required pressure balance solution is obtained from (A-8) by first making a best 
guess for the new value of K. This is done by using the "new" values of the constants of the 
motion along with the cell positions computed at the previous time step. The K selected in this 


way will clearly depend on position, and we refer to the quantity as K,,;. The true value of the 


cell boundary is, say, r° + A" where r, is our first guess or lower order iteration. If 


HUI, LIEWER, AND BOOK 


A‘) << r then an expansion r, may be performed in (A-8), where only terms linear in the 
A’s need be retained. The following recursion relation is then obtained: 

SD = LK — K2 -rva Oa, aM, --- 0/44. (A-9) 
where \.° is composed of terms linear in the various A’s and the 4’s are obtained straightfor- 
wardly by comparison with the expanded form of Eq. (A-8). If K is given, (A-9) may be solved 
for A,, i =1, 2, *++, successively because by symmetry we know that there can be no change 
in the positions of the first cell boundary , and 4; =0. (Here we adopt the convention of 
beginning the numbering of the cell boundaries at the cylindrical axis.) Of course, this values of 
K will not give a proper solution unless the position of the last cell boundary, the physical wall, 
agrees with the value specified at that instant of time. For a stationary wall the corresponding 
A,, = 0, while for a moving wall A, = 6, # 0. Hence we now construct two solutions of Eq. 
(A-9) corresponding to two values of K: 

A{} = [K@ — KO -rv~ Ma, af}, --- V/A 9, (A-10a) 


A} = (Ko) — KO = VL(A (4), A &), ae + ))/A,9}. (A-10b) 


Because Eq. (A-9) is linear in K and A,, the actual perturbed wall position can be formally 


expressed in terms of a linear combination of the solutions (A-10): 


B,, on c@aA (a) Es cf) A Cb) (A-11) 
where 
By ~ As” 
(a): oe ie cae ° 
c!?) = } AO yw’ (A-12a) 
»- A 
c) BAe (A-12b) 


Ao) a 4) 


This determination of c‘”) and c‘°) now allows us to estimate the actual cell boundary positions, 


as well as K: 
AM) _ ofA {2} t. ca ’ (A-13) 
K = cf K@ + cK), (A-14) 
20 
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The above expressions would be exact if the linearization used in the derivation of Eq. (9) were 


strictly valid and if the K;°} for all i were equal. 


The values of A,;,,; and K as determined from the above equations are treated as iterated 
quantities. Note that the iterated quantities do not depend on the choice K) and K“). At this 
stage the best guess for the cell boundary positions is r;‘!) = r, + A,,,, where A,,, is given by 
Eq. A(13). The iterated values at r;“") are use to find K({, {, and 4,42, and the recursive 
equation for the second iteration is 

A® =(K — KY -r~ Ma, a@® ---)1/4,Q). 
The procedure explained above is now repeated until the desired degree of convergence is 


obtained. 
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Appendix B 


LISTING OF THE CODE, QUEST 


The subroutines appear in alphabetical order. Member AAAI is a macro which sets the 


® 
number of grid points and species. 
t 
22 
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Appendix C 
Initial and Final Conditions for Test Problems 


ADIABATIC COMPRESSION 
Initial Conditions 


&INPUT 
GAM = 1.6666670, DRG = 0.25000000, ETA = 1.0000000, RO = 30.000000, ALFA = 25.000000, 
DENO = 0.00000000, , BZO = 2400.0000, TFACT = .0000000, BI = 100.0000, . RRIGHT = 7,0000000, 
RMIN= 1.0000000, TMID = 0.59999991E-03, TO = 10.0000000, DT = 0.59999993E-05, RWALL = 7.0000000, 
NPRINT = 25, NSTOP = 200, BTO = -2000.0000, TLBASE = 0.00000000, TRBASE = 0.00000000, 
NFRAC = 0.10000000E-03, TFRAC = 0.80000001, NPLOT = 5, 
&END, 
&KNOB, 
RJP = 2.0000000, DRJZ = 1.0000000, DRJT = 1.5000000, 
LPRINT = T, LDSP = F, LWALL = T, LPLOT = T. LBOOK = T, 
LTURB = F, LDUMP = F, 
&END, 
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SEEBIE EXPERIMENT 


Initial Conditions | 
&INPUT 
GAM = | 6666670 DRG = 0.28000000 ETA = 1.0000000, RO = 30.000000,. ALFA = 25.000000, 
DENO = -0.00000000 BZO 2400. 0000 TFACT = §.0000000, BIL = 100.00000, RRIGHT = 7.0000000, 
RMIN = 10000000 TMID = 0.99999990E-04, TO = 10.0000000, DT = 0.99999973E-08, RWALL = 7.0000000, y 
NPRINT = 10 NSTOP = 400, BTO = -2000.0000.  TLBASE = 0.00000000, TRBASE = 0.00000000, 
NFRAC = 0 |0000000E-03, TFRAC = 0.80000001 NPLOT = 100, RIP = 2.0000000, DRJZ = 1.0000000, 
DRIT = 15000000 j 
&END 
&KNOB 
LPRINT = 1 LDSP = 1 LWALL = F LPLOT = T, LBOOK = T. | 
LTURB = F LDUMP = I 
&END 
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Fig. 1 — Initial profiles of (a) magnetic field components B, and By and (b) 
plasma and magnetic pressure, taken from Ref. 7. All quantities are normalized 
to their peak values for display purposes. 
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Fig. 2 —- As in Fig. 1, after 0.3 msec of compression at 104 cm/sec, so that the 
outer boundary is at R).), = 3.9 cm. All quantities are again normalized to 
their peak values. 
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3 — As in Figs. 1 and 2, after 0.6 msec of compression, with R,,,.,, equal 
to about 15% of its initial value. 
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HUI, LIEWER, AND BOOK 
TIME (usec) 
{a) 
PROBE SIGNAL (B,) 
J 
TIME (usec) 
(b) 
Fig. 4 — Diamagnetic signal and Bg probe signal vs. time, as observed experi- 
mentally (solid trace) and calculated using QUEST, starting with the initial data 
of Fig. 1. 
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